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The proposed PINGU project (Precision IceCube Next Generation Upgrade) is expected to collect 
0(10®) atmospheric muon and electron neutrino in a few years of exposure, and to probe the 
neutrino mass hierarchy through its imprint on the event spectra in energy and direction. In 
the presence of nonnegligible and partly unknown shape systematics, the analysis of high-statistics 
spectral variations will face subtle challenges that are largely unprecedented in neutrino physics. We 
discuss these issues both on general grounds and in the currently envisaged PINGU configuration, 
where we find that possible shape uncertainties at the (few) percent level can noticeably affect the 
sensitivity to the hierarchy. We also discuss the interplay between the mixing angle 623 and the 
PINGU sensitivity to the hierarchy. Our results suggest that more refined estimates of spectral 
uncertainties are needed in next-generation, large-volume atmospheric neutrino experiments. 


I. INTRODUCTION 


The discovery of atmospheric neutrino oscillations in the Super-Kamiokande (SK) experiment in 1998 [T] marked 
the birth of what has now become the standard Zv mass-mixing scenario, where the neutrino flavor states {Vf., u^) 

are mixed with different states (j^i, 1^3) having definite masses (mi,1712,1713), via three mixing angles (012,6113,023) 

and a possible CP-violating phase 5 [ 5 ]. Several oscillation experiments have allowed to measure the three angles 9ij 
(but not yet the phase 6), as well as two independent squared mass differences among Amf^ = [ 5 ], which 

can be conventionally chosen as Sm? = > 0 and [ 3 J H] 


= ± 


Am§i -|- Am§2 
' 2 


( 1 ) 


where the sign of Am^ distinguishes the so-called normal (-I-) and inverted (—) neutrino mass hierarchy (NH and IH, 
respectively). Recent oscillation analyses with updated bounds on 0 ^- and Amfj are reported in [SHZ]. 

Global data analyses have already shown a slight sensitivity to the hierarchy through various subsets of data, 
including reactor events [ 3 ], solar events [1], atmospheric events [S] and, more recently, long-baseline accelerator 
events mi- In particular, in the observable zenith distributions of atmospheric neutrinos, sub-horizon events carry 
information on sign(ATO^) via its interference with the effective neutrino potential V in the Earth matter [ 9 ], 

V = ±V2GFNe, ( 2 ) 

where the upper (lower) sign refers to v iv), Gp is the Fermi constant, and Nf, is the electron number density. Even 
if V and V are not separated, the atmospheric neutrino event rates are not symmetric under v exchange [ 5 ], and 
thus hierarchy effects are not canceled. However, within the available data, such subleading effects have not emerged 
yet, being largely smeared out by the relatively coarse resolution in neutrino direction and energy [ 3 ]. Indeed, in 
global three-neutrino fits to atmospheric data only, the difference between the two hierarchies amounted to a mere 
Ay2 ^ Q 3 pre-SK era [ 5 ], and it is still as small as Ay^ ~ 0.9 in the latest SK Collaboration analysis m- 

However, hierarchy effects may well emerge in next-generation, large-volume atmospheric neutrino experiments, 
provided that the event statistics and the resolutions in energy and angle are high enough, as recently emphasized 
in m- In this context, the proposed ice-Cherenkov detector PINGU (Precision IceCube Next Generation Upgrade) 
[T^ . which is expected to collect O(IO^) neutrino events in a few years, appears to be a very promising project, 
and is currently under extensive investigation. Comparable goals might be reached in the proposed water-Cherenkov 
detector ORCA (Oscillation Research with Cosmics in the Abyss) [T3] and Hyper-Kamiokande [HI, while the ICAL- 
INO project (Iron CALorimeter at the India-based Neutrino Observatory) [H] aims at increasing the sensitivity to 
the hierarchy with a different approach [v and V event separation). We refer the reader to |Ib] for a recent survey of 
possible approaches to the hierarchy discrimination, by using atmospheric (and other) neutrino sources. 

In this paper we focus on PINGU, taken as a case study for analyzing the atmospheric neutrino sensitivity to the 
hierarchy with very high statistics. Several works have already addressed this task [T7H3T] . showing that PINGU 
can reach a sensitivity of at least a few standard deviations in a few years, especially in favorable conditions for 
large matter effects (i.e., for normal hierarchy and 023 in the second octant). While we confirm this broad picture 
with an independent analysis, we try to bring to surface several subtle issues related to the calculation of spectral 
shapes and to the estimate of their systematic uncertainties, which represent very peculiar challenges of high-statistics 
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data sets. Although such issues have already emerged in other fields, such as parton distribution fitting ^21 US] and 
precision cosmology [SHEl], they have only been touched upon in neutrino physics (see, e.g., [26]) and deserve renewed 
attention. In addition, we elucidate the interplay between the determination of the hierarchy and the measurement 
of 023 in PINGU. 

Our work is structured as follows. In Sec. II we describe our calculation of spectral event rates in PINGU and 
their binning in energy and (zenith) angle. In Sec. Ill we describe some general features of the statistical analysis of 
PINGU prospective data, including oscillation parameter uncertainties and normalization systematics, as well as other 
possible (correlated and uncorrelated) shape systematics, which play a relevant role in the limit of high statistics. In 
Sec. IV we analyze in detail the impact of shape systematics on the PINGU sensitivity to the hierarchy. Finally, in 
Sec. V we discuss the interplay between the hierarchy and 023 constraints. We conclude with a brief summary and 
perspectives for further work in Sec. VI. 


II. CALCULATION OF ENERGY-ANGLE SPEGTRA IN PINGU 

Our calculation of energy-angle spectra is described below. The reader not interested in details may jump to the 
last two subsections, where we discuss features of the spectra which are relevant for the subsequent statistical analysis. 


A. Notation 


The following notation is used hereafter: 


E', 9' 

E, 0 

rUE,E') 

rf{9,E) 

T 

d2$“/(dcos0'dU') 

^^ciE') 

Pc.p{0',E') 

[015 0i-|-l] 

[Ej, Ej+i] 


= flavor index (p,, e) of iz or P, 

= number of I'a + Pa events, 

= true neutrino energy and zenith angle, 

= reconstructed neutrino energy and zenith angle, 
= energy resolution function, 

= angular resolution function, 

= detector livetime, 

= effective detector mass at energy E', 

= double differential neutrino flux ($ for V), 

= ratio of double differential neutrino fluxes, 

= neutrino charged-current cross section (ct for p), 
= oscillation probability of —>■ t'/S (P for P), 

= range of i-th angular bin, 

= range of j-th energy bin. 


Note that 9' jn = 1 


and 0.5 correspond to vertically upgoing and horizontal neutrino directions, respectively. 


B. Neutrino production and propagation 

Goncerning the unoscillated atmospheric neutrino fluxes, we assume the azimuth-averaged values of 
/ {d cos 9'dE') calculated at South Pole in [28]. These fluxes are modulated by the oscillation probabilities 
Pap during neutrino propagation from the source (assumed to be a layer at = 15 km in the atmosphere) to the 
detector. For sub-horizon trajectories, matter effects are calculated up to a second-order Magnus expansion in each 
Earth shell as in [21] , along an accurate model of the electron density profile m- 

The probabilities depend, in general, on all the oscillation parameters {5m?, ±Ato^, 0i2 , 0 i3 , 023, 5). Whenever we 
need to fix “true” oscillation parameters to calculate an input spectrum for subsequent fits, we assume the following 
representative input values [5]: 


1^^ Itrue 

= 2.40 X 10"^ eV^ , 

(3) 

Itrue 

= 7.54 X 10"® eV^ , 

(4) 

sin ^13 Itrue 

= 0.0237 , 

(5) 

sin^ ^12 Itrue 

= 0.308 , 

(6) 

Itrue 

= 37r/2 . 

(7) 
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The parameter 023 is treated differently since, as discussed later, it induces large variations in the PINGU sensitivity 
to the hierarchy. Unless stated otherwise, we assume by default that it can take any true value in the range 

Sin2 023 |true G [0.4, 0.6] , (8) 

which covers both octants of 023 (as currently allowed at ~ 2 (t level in 0). Reconstructed values of sin^ 023 are 
allowed to extend beyond this range in fits. The effects of prior ranges different from Eq. Q are discussed in Sec. V. 

C. Neutrino detection 


Concerning the PINGU detector, we basically assume the preliminary characterization reported in m- We approx¬ 
imate the effective detector mass pVJ^{E') (i.e., the ice density times the effective volume for a = /r, e) by interpolating 
the E > 1 GeV histograms in Fig. 6 of [H] with the following (smooth and monotonic) empirical functions, 

pV^%{E') = 3.33 (l - , 

pV^s{E') = 3.44 (l - , 

where [pV),^] = MTon, [E'] = GeV, and the effective threshold has been set at 

EV = 1 GeV . (11) 


(9) 

( 10 ) 


The ratio of pUgff proton mass rup provides the effective number of target nucleons. 

Total charged current (CC) cross sections and for E' > E[t^^ are extracted from Fig. 14 of [12] . They are 
the sum of three contributions: (i) deep inelastic, (ii) quasi-elastic and {in) resonant scattering, where the first one 
dominates for E' above a few GeV. For simplicity, we assume that the total CC cross section for v,, is identical to the 
one for Vp at any E' > E{^^ (and similarly for antineutrinos). 


ahc{E')=a^cc{E') = acc{E') . 


( 12 ) 


The resolution functions are extracted from the 2-dimensional histograms in Fig. 14 of [T2| in digitized form Ell¬ 
in particular, in each x-axis bin having median true energy E' therein, we fit the histogram contents with gaussian 
functions having widths (Je{E') and ag{E'), 


rUE,E’) 

r^{0,0') 


1 

1 

fE-E'^ 

2 

^a%{E') 

2 

\^%{E')j 


1 

' i 

f 0-0' \ 

2 ‘ 

. - GXP 

V^a^{E') 

2 

\<^g{E')) 



(13) 

(14) 


The resulting collection of widths a^{E') and ag{E') in each bin of E' are finally fitted with the following (smooth 
and monotonic) empirical functions: 


a>^/E' = 0.266/(F'°'i^i - 0.604) , (15) 

ctI/F' = 0.369/(F'°-2^^ - 0.508) , (16) 

cr^ = 3.65/(F'1°5 -h 5.00) , (17) 

= 1.88/(F'° ®23 -p 1.93) , (18) 


where [E'] = GeV and \ag\ = rad. These approximations capture the main features of PINGU as described in [12] . 

Figure 1 shows the ±lcr resolution bands for Vp in the plane charted by log2o(7?VGeV) versus 0'/tt (left plot) or 
versus cos0' (right), in the intervals E' G [1, 40] GeV and 0 G [0.5, 1].^ The horizontal bands correspond to Ea^iE') 
for E' = 3, 10, and 30 GeV, while the three curved, vertical bands correspond to E(jg{E') for 0' /tt = 0.6, 0.75 
and 0.9. The resolution functions largely smear any spectral feature when passing from true to reconstructed variables, 
{E\ 0') —)■ (F, 0), the more the lower the energy. 


^ As usual, upgoing events {O/tt ~ 1) correspond to the left of the zenith scale, and horizontal events {O/tt ~ 0.5) to the right. 
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FIG. 1: Widths (at ±1 (t) of the resolution functions in energy and angle for r/^-like events in PINGU, in terms of logarithmic 
energy versus the zenith angle (left plot) or its cosine (right plot). See the text for details. 


Figure 1 also illustrates three advantages of using the zenith angle rather than its cosine. The first is that the 
angular resolution bands (which provide a rough idea of the appropriate zenithal binning) are obviously symmetric in 
9' (left plot) but not in cosd' (right plot), where they are squeezed towards the upgoing directions {cos9' —>• —1). The 
second is that, compared with the full sub-horizon range ff'/rr G [0.5,1], the interesting angular fraction subtending 
the dense Earth core (9'/7r G [0.816,1]) is as large as 36.8%, while it would be squeezed by a factor of about two 
(16.2%) in terms of cos0'. The third is that, by using cos 9', one would expand the nearly horizontal part of the 
zenith spectrum (cos 9' > —0.5) which, despite being weighted by higher atmospheric fluxes [3H], is less interesting 
for hierarchy discrimination, due to smaller matter effects at shallow depth in the Earth’s mantle. 

Concerning the energy, we remind that in normal (inverted) hierarchy, matter effects for neutrinos (antineutrinos) 
are particularly enhanced around E' ~ 2.5-3 GeV and E' ~ 6-10 GeV (corresponding to resonance effects in the 
core and in the mantle, respectively), as well as for intermediate energies where mantle-core interference effects occur 
(see [HI] and refs, therein). Although the low-energy range E' ~ 0(1 — 10) GeV contains most of the hierarchy 
“signal,” it is useful to extend the analysis to few tens of GeV (or more), for at least two reasons: (1) the high-energy 
spectrum is better experimentally resolved and is largely hierarchy-independent, so it can help to “fix” some floating 
parameters in the fits; (2) due to the relatively poor energy resolution at low energy, hierarchy effects may “migrate” 
well above ~ 10 GeV in reconstructed energy. In any case, to avoid “squeezing” the most relevant low-energy range, 
it is useful to adopt a logarithmic energy scale, as in Fig. 1. Summarizing, we shall use the zenith angle (instead of 
its cosine) and a logarithmic energy scale for representing PINGU event spectra. 

D. Unbinned event spectra 


In atmospheric neutrino experiments, the detected neutrino events are usually organized in terms of energy and 
zenith angle (or related variables). The double differential spectra of iV“ events induced in PINGU by both and 
pQ, as a function of the true energy E' and zenith angle 9', can be cast in the form 

dcos 9'dE' 


Too dcos 9'dE' 


P^(9',E') 


(19) 


where the prefactor in square brackets does not depend on the oscillation parameters, while the last factor P°‘ is a 
linear combination of the relevant oscillation probabilities: 


pa _ 


$/5 


-b 


crcCo , ^ 

^ n/CV \ ^ 0Ct 


4>“ CTCC 




( 20 ) 


where jd ^ a and the first (second) term in brackets is due to v (ly), respectively. In Eq. (19) we have assumed a 
priori azimuthal averaging, hence the 27r factor; see m for an approach including azimuth dependence. This issue 
and other related approximations will be discussed in Sec. II G. 

The spectra in terms of reconstructed variables (E, 9) are obtained by convolving the ones in Eq. (19) with the 
resolution functions. 


d9dE 


n2-K 


sin9'd9'r'^(9,9') 


\{E,E') 


d'^N° 


d cos 9'dE' 


(21) 


where the change of variable cos 9' —)• 9' has been applied, as discussed at the end of the previous subsection. 
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E. Binned event spectra 


The number of events in the rj-th bin is obtained by integrating the r.h.s. of Eq. (l2lj) over the bin area 
[6i, 9i+i\ 0 [Ej, Ej+i], By changing the integration order (see [311 [32]), the resulting quadruple integration can be 
reduced to a double one: 

( 22 ) 

(23) 

where the functions 'w^{x) are defined, for {n,x') = {i,0') and {j,E'), as: 


r^i+i 


rEi. 


r2-K 


7V“ = 




d9 dE sm9'd9'r^{9,9') dE'TE{E,E')- 

j E- Jo JE' dcos9 dE 


nOO 

sin6»'d6»' / dE'wf{9')w'J{E') 

Je' 


thr 

r2 Aja 


d^N° 


d cos 9'dE' ’ 


/ /\ t f ( ®ra+l 2: 

''''nix ) = ^ erf ' 


with erf(a::) defined as in 


V2o 


2 V ^2(7“ 


(24) 


erf(a:) = 

'Je 


dte 


(25) 


In the limit of perfect resolution (tr" —>■ 0), the curve Wnix') becomes a top-hat function in the interval [xn, a;n+i]- 
For finite resolution, the “top-hat” shape is smeared and extends beyond this interval for a few cr“’s. However, for 
numerical purposes, Wn{x') practically vanishes beyond — 4cr“, Xn+i + 4(t“], so that the double integral domain 
in Eq. (23) can be just taken as the ij-th bin range “augmented” by ±4crg and ±4 (t^. 

As previously argued, we actually adopt a logarithmic energy variable. 


A = logio(E7GeV) , 


(26) 


(and similarly for the reconstructed energy E), so that 


^277 




sin6»'d6»' / dA£;'ln(10) w“(6»'X(E') 


d^N° 


d cos 9'dE' 


(27) 


Finally, we consider reconstructed energies in the interval E € [1, 40] GeV, namely, in the logarithmic range 
log 2 o(E/GeV) S [0, 1.6], that we divide into 16 bins. We also divide the range of sub-horizon reconstructed an¬ 
gle, S/tt = [0.5, 1], into 10 bins. With this choice, the bin widths are smaller than the typical resolution widths in 
Fig. 1 (so as to avoid additional smearing from binning), but large enough to contain a signihcant number of events 
after a few years of exposure (so as to apply Gaussian, rather than Poissonian, statistics). The calculation of AT. via 
Eq. is performed through Gauss quadrature routines, which have been checked to yield numerically stable results 
up to the third significant hgure, even in bins where the integrand oscillates rapidly via P“. 


F. Qualitative discussion of spectral shapes 


Figure 2 shows the main ingredients of typical PINGU spectra calculations (first three couples of panels from the 
left) and final spectra (last couple of panels on the right), where the upper and lower panels refer to muon events 
(a = /i) and to electron events (a = e), respectively. The adopted ranges and bins have been discussed in the previous 
subsection. For definiteness, we have assumed normal hierarchy (NH), sin^ 023 = 0.5, and the remaining oscillation 
parameters as in Eqs. in any case, the graphical results would appear qualitatively similar for different choices. 

The units and the color scale are arbitrary: in each panel, the darkest color corresponds to the bin with maximum 
contents, while lighter shades refer to lower contents, down to total white for almost empty bins.^ 


The leftmost panels in Fig. 2 show the product namely, the oscillation-independent prefactor in Eq. (191. 

This prefactor is suppressed at high energy by $“crg(;; ^ and at low energy by the small value of VJg, with 


^ Concerning the absolute event rates, for the specific oscillation parameters chosen in Fig. 2, we estimate a total of 1.9 X 10^ muon and 
1.4 X 10'^ electron events per year in PINGU. The total statistics can thus reach 0(10®) events in a few years, as already noted. 
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V“, 0“ ogc Pnh (^23=0.5) (unsmeared) (smeared) 



0/71 0/71 0/71 0/71 


FIG. 2: Breakdown of the main ingredients of typical muon (upper panels) and electron (lower panels) event spectra in energy 
and angle in PINGU. From left to right: oscillation-independent factor oscillation-dependent factor P“, their 

product (unsmeared event spectrum), and final results (smeared event spectrum). See the text for details. 


a maximum in the few GeV range, which is interesting for matter effects. Unfortunately, the atmospheric neutrino 
flux peaks at the right energy but in the wrong direction, i.e., at the horizon (S/tt —)■ 0.5), where matter effects 
vanish. In this sense, atmospheric neutrinos are not “optimal” for seeking hierarchy effects, which occur mainly in a 
tail (rather than at the peak) of the event spectrum. 

The next couple of panels in Fig. 2 show the oscillation-dependent factor in Eq.(20). In the upper panel, the 
factor shows large variations, whose shape is reminiscent of the oscillograms related to disappearance [33]. In 
particular, a large disappearance “valley” (the first oscillation minimum) extends from the upper left corner to the 
lower right margin of the panel. Conversely, in the lower panel, the factor P® shows much milder variations, since the 
iZg disappearance and appearance probabilities are largely suppressed by the smallness of sin^ ^ 13 . 

The third conple of panels shows the binned product of the factors and P“, which is proportional to the 

unsmeared spectrum of events in Eq. (19). An oscillatory structure is still visible in the central part of each panel, i.e., 
for slanted trajectories and for E ^ few GeV. These structures, however, are largely suppressed around the vertical 
upgoing direction, where the atmospheric flux is lower. 

Einally, the rightmost couple of panels shows the observable, smeared spectra of ^ and e events, including resolution 
effects as in Eq. ( |27| . The oscillating structures appear to be largely suppressed, except for remnants of the large 
disappearance valley, carrying the dominant information about the oscillation parameters (|Am^|, sin^d 23 )- The 
smeared spectra in inverted hierarchy (not shown) would be visually indistinguishable from the ones in normal 
hierarchy in Eig. 2. Indeed, as discussed in the next Section, hierarchy effects emerge as relatively smooth and 
subdominant modulations, at the level of a few percent, in the left tail of the zenith spectra and for intermediate 
energies. It is thus imperative to assess how accurately we know the shapes of the neutrino event spectra in PINGU 
and, in general, in large-volume atmospheric experiments. 
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G. Qualitative discussion of shape uncertainties 


In a few years, PINGU is expected to collect O(IO^) events which, distributed in O(IO^) bin contents N^, can 
constrain the energy-angle spectral shape with a statistical accuracy 1 /^ few%. However, the shape may also 
be affected by comparable systematic errors stemming from the oscillation-independent factors (I4g$“crg(;;), 

the oscillation-dependent ones (P“), the resolution functions (r|; g), and the event integration into bins.^ 

The effective detector volume V^, which carries information about the event detection and reconstruction efficiency, 
has been estimated in PINGU via Monte Garlo (MG) simulations [T^]. In particular. Fig. 6 therein shows its E' 
dependence, with fluctuations suggestive of finite MG statistics. In general, the event collection efficiency depends 
also on 9' (due to the cylindrical configuration) |II5j and, to some extent, on the azimuth angle (j)' (due to unavoidable 
anisotropies of the detector and inhomogeneities of the ice). Apart from normalization and finite MG statistics errors, 
it is reasonable to think that the V),g’s are affected by shape uncertainties SV^{E', 9', (j)') at the (few) percent level, 
possibly different for p and e, especially in the relevant range E ^ few GeV, where the lUff'® change rapidly.^ 

The atmospheric neutrino fluxes $ also depend, in general, on the three variables {E', 9', </)'). Overall normalization 
errors are large [O(I0)%], but are also highly correlated between {p, e) or {ly, p), and can thus be reduced to a (few) 
percent by taking appropriate ratios. At this level of accuracy, a number of irreducible shape uncertainties are 
also known to be present. First-order shape variations have often been parametrized via “tilts” of both the energy 
spectrum index and the zenith-angle distributions (see, e.g., |371188) 1. This simplification is generally adequate when 
the investigated signal is sizable, as it happens for atmospheric disappearance due to the dominant (|Am^|, ^ 23 ) 
parameters, but may be too coarse for subdominant signals. A dedicated Workshop recognized, a decade ago, 
the need for a more refined characterization of atmospheric flux (and other) uncertainties, in view of future high- 
statistics experiments seeking subleading effects [26]. As a relevant follow-up, the authors of |39| broke down the main 
atmospheric flux uncertainties into Ns = 26 independent error sources within a one-dimensional propagation model 
for the Kamioka site, and studied their associated effects at the (few) percent level on the energy-angle spectra.® It 
would be desirable to repeat this relevant study in a three-dimensional flux model for the South Pole site, possibly 
within two or more independent simulations, in order to identify an analogous set of Ns systematics and associated 
flux variation functions Ss^°‘{E',9 'to be included in data fits. At present, it is legitimate to assume that the 
shapes of the $“(i?', 9', ^') functions are not known to better than the (few) percent level. 

Concerning the calculation of the oscillation probabilities Pap, the main uncertainties are related to variations of 
the mass-mixing neutrino parameters, particularly of the dominant ones (|Am^|, ^ 23 ); such variations are known to 
reduce the sensitivity to hierarchy effects m- To a much lesser extent, the Pap are affected (in the atmosphere) by 
uncertainties in the production height distribution, and (in matter) by uncertainties in the electron density profile. 

The CC cross sections are poorly known in the (few) GeV range, especially when deep inelastic scattering is 
not dominant. Uncertainties on total and differential cross sections affect the spectral normalization and shape, 
respectively, at a typical level of few %. The shape uncertainties may affect high-statistics estimates of the dominant 
parameters (Am ^,023 ) m- A fortiori, such uncertainties cannot be ignored in the extraction of subdominant effects. 

Detector response uncertainties add upon cross section errors in the determination of the energy-angle resolution 
functions, which characterize the probability distributions of the observable (reconstructed) kinematical parameters 
around the unobservable (true) ones. The resolution functions are reported in |12] (see Fig. 8 therein) in terms of 
E\ but they may also depend on (0', (j)') via detector response anisotropies. In the absence of a calibration neutrino 
beam, the resolution functions can only be simulated, and their centroids and shapes are unavoidably affected by 
both cross section and reconstruction uncertainties. 

Finally, the multi-dimensional integration into binned spectra may be a source of numerical errors in itself. For 
instance, in order to keep the calculations manageable, we have reduced the number of nested integrations, by averaging 
out a priori the azimuth angle, the production height, and the energy and angle(s) of the outgoing leptons (Sec. II D). 
We have also assumed gaussian resolution functions in order to apply, in each bin, the reduction in Eq. (231. Of course. 


these approximations can be avoided by a brute-force integration over all the relevant variables, including the (true 
and reconstructed) neutrino and lepton energies and directions. This approach (which is well beyond the scope of this 
paper) becomes impractical, however, when one must repeat the calculations by varying many systematics in real or 
prospective data fits. Note also that the oscillation probabilities may vary wildly over typical energy-angle bin ranges, 
making the calculations quite demanding in terms of integrand grid sampling and numerical accuracy. 


® Although some spectral shape errors have been previously considered in the literature 1111112111714211 . we think it useful to present a 
more extended and self-contained discussion herein. 

^ For the sake of comparison, in a laboratory neutrino experiments such as Borexino m, deviations of the effective fiducial volume from 
its quasi-spherical shape are estimated to be at the few percent level from point to point [0(10) cm over a few-meter diameter]. 

^ See also m for an independent error assessment. 




Alternatively, one might replace multi-dimensional integration into bins by a full Monte Carlo approach, randomly 
following the entire process of production, propagation, interaction, reconstruction and binning of events.® In order 
to keep the statistical MC error at subpercent level, each bin should collect no less than 0(10"^) simulated events, 
which implies 0(10®“^) events for each MC spectrum with O(IO^) bins. Moreover, the MC spectrum must be 
repeatedly calculated, by sampling the probability distributions of the floating variables in the fit (including the 
oscillation parameters and the systematics uncertainties), leading to no less than 0(10®) generated events, i.e., to a 
MC statistics several orders of magnitude higher than the real data sample of 0(10®) events. It is not obvious (at least 
to us) that this goal, and the desirable subpercent statistical MC accuracy in each bin, can be practically reached. 
We thus argue that any method of calculation of event rates can also induce (largely uncorrelated) residual errors in 
each bin, at about the percent level. These additional numerical uncertainties may not be totally negligible in the 
statistical analysis of next-generation atmospheric neutrino experiments such as PINGU. 

In conclusion, we have discussed qualitative arguments in favor of many possible sources of (few) percent uncertain¬ 
ties on the shape of the energy-angle distributions in PINGU. Each of these sources gives extra freedom to adjust the 
event spectra in data fits, thus reducing the sensitivity to spectral differences induced by hierarchy effects. Although 
each reduction may be small, their sum may become noticeable. In the following Section we shall study, in a more 
quantitative way, the progressive impact of some of the above error sources in PINGU. 


III. STATISTICAL APPROACH TO THE HIERARCHY SENSITIVITY IN PINGU 

In this Section we discuss our statistical approach to the hierarchy sensitivity in PINGU. We define the methodology 
used to deal with systematic errors, and discuss some subtle issues emerging in the limit of very high statistics. We 
then consider an increasingly large set of errors, including those coming from oscillation parameter and normalization 
uncertainties, from (known and unknown) correlated shape systematics, and from possible residual uncorrelated errors. 


A. Methodology for systematics, and discussion of high-statistics limits 


Our statistical analysis of the and e event spectra in PINGU is based on a y® approach. As argued in [32], 
a good metric for the sensitivity to the hierarchy (despite its discrete nature) is given by the marginalized Ay® 
difference between the true hierarchy (TH) and the wrong hierarchy (WH), where (TH, WH) refer to either (NH, IH) 
or (IH, NH). The parameter AU = a/A y® can thus be taken as a shorthand for the effective number of standard 
deviations separating the TH and WH hypotheses. 

The TH and WH spectral event rates are defined as 


R%{pk) = 

Pk) 

(28) 

T 

R%{pk) = 

iV“(WH; Pk) 

(29) 

T 


where T is the detector live time, the pk are the (oscillation and systematic) fixed parameters in TH, while pk are 
the corresponding floating parameters in WH.^ The “theoretical” WH hypothesis is tested against the “experimental 
data” represented by the TH event rates, which are affected by statistical errors decreasing as Vt, 


Vt 


(Ic^) > 


and by systematic errors on the parameters pk , 


Pk ± <Jk (Icr) . 


(30) 


(31) 


^ However, to our knowledge, there is no known MC code including, at the same time, both atmospheric neutrino production from cosmic 
rays and neutrino interactions in a South Pole detector, since the associated MCs have been developed by different groups of researchers. 
^ Since we neglect seasonal variations and take average fluxes from EH], the event rates are constant in time. 
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The function is defined as 


Ax = niin 

Pk 


10 16 

EE E 

i—l j — 1 oc—fi, e 


(Pfc) - i?“ 

iW+W) 


k 


Pk -Pk 
^k 


( 32 ) 


where the second term represents the sum of penalty functions for the nuisance parameters pk (assuming gaussian 
errors cife). Special cases for the penalties are: (1) Cfe = 0, which corresponds to having the fc-th parameter hxed as 
Pk = Pk', and (2) Ufe = oo, which corresponds to unconstrained values of pk- 

In the above equation, the first term includes possible uncorrelated errors which (as argued in Sec. II G) may 
stem, e.g., from finite MC statistics in each bin, or from residual systematics. There are however, deeper motivations 
to include such uncorrelated errors, as noted in [33]. In fact, the above Ay^ function entails two strong assumptions: 
(1) that we know all the possible sources of correlated systematic parameters pk', and (2) that we know exactly the 
effect of each pk variation on each binned rate . In the limit of very large statistics (s“ —> 0), these two assumptions 
would lead to Ax^ —t oo for Uij = 0, unless Rfj{pk) = Rfj{pk) for some very peculiar values otpk- However, in general, 
no combination of pk can exactly reproduce a generic data set . In other words, the assumptions that we do know 

all the systematics and all their spectral effects leads to the paradoxical conclusion that the sensitivity Ng- = a/A y^ 
grows indefinitely with Vt, without reaching any reasonable, systematics-limited plateau [43] . 

This situation is basically unprecedented in neutrino experiments, usually characterized by relatively low statistics.® 
To solve the paradox, one should admit that the knowledge of spectral systematics may be either incomplete or 
inaccurate to some extent, and try to deal with the residual ignorance. One possibility is to render the spectra more 
flexible, by including additional families of admissible spectral deviations via extra parameters pk, constrained by 
dedicated considerations or educated guesses. This method has been explored, e.g., in the context of fits to parton 
distribution functions [22l|23| and to precision cosmological data [24l|25|.® Another possibility is to parametrize our 
ignorance by allowing additional uncorrelated errors Uij of reasonable size in each bin, which lead to hnite Ay^ values 
in the limit of infinite statistics |33|. In our analysis, we shall study in sequence the effects of both approaches (which, 
in a sense, try to deal with “errors on errors” [46]1. 

The x^ minimization procedure is also nontrivial for large numbers of bins i x j and of floating parameters pk- 
In the context of PINGU, a brute-force scan of the {pk} parameter space through a fixed sampling grid would be 
prohibitively time consuming, also because the hierarchy sensitivity can be as large as Ncr ^ 0(10) in favorable cases, 
implying that the pk distribution tails must be sampled at the same level. Alternatively, one might explore the large 
{Pk} parameter space via a Markov Ghain Montecarlo (MCMC) method |3Z|. This method (as the brute-force scan) 
does not make any hypothesis on the functional form of the functions Rfj{pk), and can work also for generic (non¬ 
quadratic) penalty functions for the pk- However, we have found that, at least in our implementation of the MCMC 
[32l|47] for the PINGU analysis, the cases with high sensitivity ^ 0(10) are numerically too demanding and time 
consuming; therefore, we have abandoned this approach. We have finally chosen to use a method mainly based on the 
“pull approach” [48] . which increases enormously the minimization speed and stability, at the price of of assuming a 
first-order expansion of the Rfj{pk) under small deviations of (some) parameters pk around the “true” values pk- 


Pk =Pk+^k 


RUPk) ^ RtApk) + 


dRUPk) 

dpk 


(33) 


Pk—Pk 


In the context of our PINGU analysis, the linearization of the rates is actually exact for some parameters (e.g., 
normalization errors), and is a reasonably good approximation for all the other parameters, with the only exception 
of the mixing angle 023 and the CP-violating phase S, whose dependence cannot be linearized at all. For this reason, 
for any given choice of TH parameters (sin^ 023, ^), we do scan the WH parameters (sin^ 023, ^) over a grid sampling 
the full range [0, I] G [0, 27r]. For each (sin^023,5) point of the the grid, we numerically calculate the derivatives in 
Eq. (33) by taking finite differences at zL2ak, then minimize the x^ analytically over the linear(ized) pk variations 


and finally find the absolute minimum by scanning the whole (sin^ 023,(5) grid. 


® Actually, current short-baseline reactor experiments with huge statistics (million events) are now facing such challenges in the charac¬ 
terization of spectral shape systematics, see El and references therein. 

^ Pushing this approach further, the likelihood could eventually be minimized over a functional ensemble (via path integral techniques) 
rather than over a discrete pk ensemble, see m]. 
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B. Systematics due to oscillation and normalization uncertainties 


The most obvious sources of systematic errors are due to: (1) the calculation of oscillation probabilities and (2) 
absolute and relative normalizations. Concerning the first, we attach the following Icr fractional uncertainties to the 
central values in Eqs. © and © 0: 

a{Am^) = 2.6% , (34) 

cr(sin2 6 »i3 ) = 8.5% . (35) 

The parameters and sin^ di 2 are kept fixed as in Eqs. © and ©, respectively, since their errors induce negligible 
effects in the PINGU analysis. As already discussed, the true parameter S is fixed at 3n/2 as in Eq. ©, while for the 
wrong hierarchy it is left free in the range [0, 27 r]. The true parameter sin^ 023 is chosen in the range [0.4, 0.6], while 
for the wrong hierarchy it is left free in the range [0, 1].Finally, we add a reasonable 3% error on the electron density 
in the Earth’s core, 


cr{Nf.) = 3% (core) , 


(36) 


to account for uncertainties in its chemical composition.^*^ 

Concerning the absolute normalization, we attach an overall 15% error /jy to all the event rates, accounting for 
fiducial volume, flux and cross section uncertainties, 

^ (1 + M , aifN) = 0.15 . (37) 

The relative normalizations between the /r and e rates, and between the v and V components of the rates, are allowed 
to differ, respectively, up to 8% and 6% at Icr (which represent values in typical ranges |THI|38])^^ 


Uf, 


+ hfR) \ 
RUl - kfR) ) 


aifn) = 0.08 , 


(38) 




(t(/,) = 0.06 . 


(39) 


Note that, at this stage, we are not consider further systematics, either correlated {pk) or uncorrelated (m“ ), which 
could give more “flexibility” to the spectral shape. 

Let us discuss two examples of best-fit spectra in a “favorable” and “unfavorable” case for PINGU, including the 
above uncertainties. In general, the case of true normal hierarchy is more favorable, since the corresponding matter 
effects are stronger for neutrinos, which are enhanced by a larger cross section than antineutrinos. Cases with sin^ 023 
in the second octant are also more favorable, since this parameter modulates the amplitude of the relevant oscillation 
channel P^e, which embeds large matter effects. Conversely, cases with true inverted hierarchy and sin^ 023 in the 
first octant are typically less favorable for hierarchy discrimination. 

Figure 3 shows the absolute spectra R^j{pk) and R^jipk) in terms of events per bin (upper and lower left panels, 
respectively) for the favorable case of true normal hierarchy and sin^ 623 = 0.6. The middle panels show the statistical 
differences between such spectra and the corresponding best-fit spectra in inverted hierarchy, R'ljipk) and Rfj{pk), 
after marginalization over the pk systematic parameters (oscillation and normalization uncertainties) described above. 
As pointed out in several papers [iiiiiaiiMi], statistical differences, up to l-2tT in some bins, appear in the energy- 
angle region where matter effects are generically large, and even beyond (due to smearing); the differences typically 
change sign by changing flavor and, for a given flavor, they also change sign in the energy-angle plane. Such patterns 
of statistical deviations should thus provide useful cross-checks, provided that they are not spoiled by systematic 
effects. The right panels show the same differences, expressed in terms of percent deviations, reaching a few % for 
muon event spectra and twice as much for electron event spectra. Hypothetical systematic shape deviations at the 
5-10% level, “equal and opposite” to those shown in the right panels of Fig. 3, would basically cancel the hierarchy 
difference, strongly reducing the related PINGU sensitivity. A nonegligible reduction can still be expected, however, 
for smaller shape deviations at the (few) percent level. 


The typical difference between the Ne values in the mantle and in the core is ~ 6% m- 

In principle, due to oscillations, one should distinguish between a relative flux error at the source, and a relative ^ij/^ij niis- 

identification error at the detector. However, we have verified that these errors are highly correlated in the fit results, and their merging 
in a single error source is a justified approximation within the present work. 
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FIG. 3: Case of true normal hierarchy and sin^ O 23 = 0.6. Left panels: absolute event spectra (top: events; bottom: e events). 
Middle panels: statistical deviations with respects to the best-ht spectrum in the wrong (inverted) hierarchy, marginalized over 
oscillation and normalization systematics only. Right panels: the same deviations in percent values. 



FIG. 4: As in Fig. 3, but for true inverted hierarchy and sin^ 623 = 0.4. 
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Figure 4 is analogous to Fig. 3, but is obtained for the “unfavorable” case of true inverted hierarchy and for 
sin^ 023 = 0.4. In comparison with Fig. 3, the results in Fig. 4 show the following features: (1) the left panels 
are basically indistinguishable, confirming that the hierarchy discrimination can only emerge from careful spectral 
analyses and not “by eye;” (2) the deviations in the middle and right panels are, as expected, generally smaller (and 
with opposite sign) with respect to those in Fig. 3, making systematic shape deviations at the few % level more 
dangerous for hierarchy discrimination; and (3) the detailed patterns of deviations in Fig. 4 are somewhat different 
from Fig. 3, especially for fj, spectra, where the deviations can extend to relatively high energies and change sign twice; 
these features may be traced to the fact that the wrong (NH) /i spectrum tries to fit the true (IH) spectrum also via 
deviations of the dominant oscillation parameters, which induce a mismatch in the region of the first minimum for 
disappearance. 

In conclusion, Figs. 3 and 4 show, once more, the necessity to investigate the impact of shape systematics at the 
(few) percent level in a wide energy-angle region, not necessarily restricted to nearly upgoing events at few GeV. For 
this purpose, following the discussion in Sec. II G, we shall first include “known” systematics directly affecting the 
shape, and then try to parametrize “unknown” uncertainties. 


C. Adding energy-scale and resolution width uncertainties 


Uncertainties in the double differential cross-section, as well as in the detector energy-angle reconstruction, even¬ 
tually affect the shapes of the resolution functions. For the sake of simplicity, we consider only gaussian resolution 
functions [see Eqs. (13) and 0], which may be affected by two kinds of systematics: biases in the centroid, and 
fluctuations in the width. More complicated (e.g., skewed) variants might be considered for nongaussian cases. 

Following recent PINGU presentations (see, e.g., [35]), we assume that the true energy centroid may be biased by 
up to 5% at Icr, 


E' ^ E'{1 + /e) , aifE) = 0.05 . (40) 

We actually include two such energy scale errors, for p, and e event independently, while we neglect possible directional 
biases for such events, which are expected to be smaller (and are usually undeclared in PINGU presentations). 

Concerning the resolution widths, on the basis of our histogram fitting procedure described in Sec. II C, as well 
as on fluctuations in the PINGU own evaluation of widths |35|, we estimate that they may vary up to 10% at lu, 
independently of each other: 


r“^r“(l + /“) , a(/“) = 0.1 , (41) 

where a = /r, e and z = E, 0. The allowance for slightly wider or narrower resolution functions is a relevant degree 
of freedom in the fit, given the role of smearing effects in determining the observable spectral shapes. 


D. Parametrizing residual correlated systematics with polynomials 

The previous energy scale and resolution width errors do not exhaust the (presumably long) list of shape systematics. 
As argued in Sec. II G, uncertainties in the effective volume and in the reference atmospheric fluxes may also lead to 
an entire set of (few) percent deviations as a function of energy and angle, which are not necessarily well known or 
under good control. Usually, this ignorance is parametrized in terms of first-order deviations (“tilts”) of the spectra 
but, in the context of PINGU, one should allow for further (smooth) nonlinear deviations, which may be more crucial 
for hierarchy discrimination, as shown by Figs. 3 and 4. Nonlinear systematics at (few) percent level are known to 
affect, e.g., atmospheric neutrino flux shapes in both energy and angle [35], as discussed in Sec. II G. 

In the absence of a detailed study of such residual shape systematics in PINGU, we provisionally assume that the 
observable PINGU spectra have small, additional functional uncertainties, parametrized in terms of polynomials. In 
particular, we rescale the abscissa and ordinate variables of the spectra in Figs. 1-4 as 

X = 401^-3 , (42) 

2 / = 1.251ogio(U/GeV)-l , (43) 


so that they range within 


( 2 ;, y) G [—1, +1] G [—1, +1] . 


(44) 
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We assume that the rates i?“ may be subject to generic h-degree polynomial deviations in {x, y) of the kind: 


R° 


1 + 


h 

Y ^ 

/ j ^nm 

m+n>0 




(45) 


where Xi and yj are the midpoint coordinates of the ij-th bin. The coefficients are allowed to float around a null 
central value within representative errors, that we choose as 


a(c“^) = 10-2 X 


1.5 (default) , 

3.0 (doubled errors) , 
.75 (halved errors) , 


(46) 


so as to cover cases with shape systematics at various levels (percent, few percent, subpercent). 

We shall consider polynomials with /i = 1, 2, 3 and 4, corresponding to linear, quadratic, cubic, and quartic 
deformations in x and/or y. Such polynomial shape deformations are easily implemented in the statistical analysis, 
since the rates i?“ are linear in the additional pull parameters which provide from 4 (linear case) to 28 (quartic 
case) extra degrees of freedom in the fit.^^ By construction, the ±lcr excursion of a generic term is thus 

contained within ±1.5% in the default case [and, similarly, within ±3% and ±0.75% in the other two cases of Eq. (46)]. 
Although several terms might add up to much more than ±1.5% in the default case, such a freedom is never 

really exploited in the fit, leaving the best-fit spectral deformation at a relatively small, few-percent level, as we shall 
comment in the next Section IV. 


E. Parametrizing residual uncorrelated systematics 


It is legitimate to posit (as argued in Sec. II G) that our knowledge of the systematics, no matter how detailed, may 
be incomplete, leaving residual uncorrelated errors it“ in each bin from various sources (including finite statistics effects 
in atmospheric, cross section, and reconstruction MC simulations). We shall consider three simple, representative cases 
for the residual (uncorrelated) fractional uncertainties in each bin, namely: 


^ = 10 


-2 


1.5 (default) , 

3.0 (doubled errors) , 
.75 (halved errors) . 


(47) 


This completes our list of uncertainties, which will be progressively included, via Eqs. (32) and (33), in the following 
statistical analysis. 


IV. STATISTICAL ANALYSIS OF THE HIERARCHY SENSITIVITY: RESULTS 

Figure 5 shows the PINGU sensitivity to the hierarchy, in terms of standard deviations separating the true mass 
hierarchy (top: NH; bottom: IH) from the wrong mass hierarchy, as a function of the detector live time T in years. 
The bands cover the fit results obtained by spanning the range sin^ 023|true G [0-4, 0.6]. The abscissa is scaled as VT, 
so that the bands would grow linearly in the ideal case of no systematic errors (not shown). From left to right, the 
fit includes the following systematic errors: oscillation and normalization uncertainties, energy scale and resolution 
width errors, polynomial shape systematics (with up to quartic terms), and uncorrelated systematics, as defined in 
Sec. III. The last two error sources are kept at the default level of 1.5%. With only normalization and systematic 
errors, Wr grows almost linearly in VT, i.e., the experiment is not limited by these systematics, even after 10 years of 
data taking. However, the progressive inclusion of correlated shape systematics, both “known” (resolution scale and 
widths) and “unknown” (ad hoc polynomial deviations), and eventually of uncorrelated shape systematics, provide 
a suppression of N„, whose estimated ranges increase more slowly than -\/T. The typical effect of all the systematic 
shape errors in the rightmost panels is to decrease the 5-year (10-year) PINGU sensitivity by up to ~ 35% (~ 40%), 
with respect to the leftmost panels in Fig. 5. 


Note that, in Eq. | |45| , the terms Cqq are omitted, since they correspond to the overall normalization errors of Sec. Ill B. Note also that 
the linear terms Cjq and Cq^ parametrize the usual spectral “tilts” discussed in Sec. II G. 
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FIG. 5: PINGU sensitivity to the hierarchy {Ncr), for either true NH (top panels) or true IH (bottom panels), as a function of 
the live time T in years. The abscissa is scaled as \/T, so that the sensitivity bands (which span sin^ 023|true G [0.4, 0.6]) would 
grow linearly for purely statistical errors. From left to right, the fit includes the following systematic errors: oscillation and 
normalization uncertainties, energy scale and resolution width errors, polynomial shape systematics (with up to quartic terms) 
at the 1.5% level, and uncorrelated systematics at the 1.5% level, as defined in Sec. III. 


Table I reports numerical results for the same fit of Fig. 5, with a breakdown of the polynomial shape systematics 
(from linear to quartic deviations). It can be seen that most of the sensitivity reduction due to polynomial shape vari¬ 
ations is already captured at the level of linear and quadratic parametrization, with higher-degree terms contributing 
a small fraction of la. Although each polynomial term can typically contribute to a ±1.5% deviation by 

construction, their sum ^ yields typical deviations of about ±1% (±2%) for (e) events in the fit, except 

for the case of NH in the second octant, where they can become twice as large (but where is also large). In 
conclusion, reasonable shape uncertainties at the (few) percent level may produce a noticeable overall effect on the 
PINGU sensitivity, although none of them appears to be crucial in itself. 


TABLE I: Reduction of the PINGU sensitivity to the hierarchy (expressed in terms of No- range for sin^ ^23 € [0.4, 0.6]) dne to 
the progressive inclusion of various shape systematics, for 5 and 10 years of exposnre. Gorrelated polynomial and uncorrelated 
systematic uncertainties are taken at the default level of 1.5%. See the text for details. 


5-year sensitivity No- 

10-year sensitivity No- 

Errors included in the fit 

True NH 

True IH 

True NH 

True IH 

Stat. -I- syst (osc.-l-norm.) 

4.23-12.3 

3.34-5.64 

5.82-16.1 

4.49-7.64 

+ resolution (scale, width) 

3.31-9.76 

2.95-4.37 

4.54-12.9 

4.00-5.94 

+ polynomial (linear) 

3.14-9.17 

2.86-4.16 

4.23-11.9 

3.81-5.49 

+ polynomial (quadratic) 

3.01-8.29 

2.69-3.88 

3.93-10.6 

3.47-5.05 

+ polynomial (cubic) 

2.98-8.26 

2.67- 3.84 

3.87-10.5 

3.42-4.94 

+ polynomial (quartic) 

2.95-8.12 

2.64-3.79 

3.82-10.3 

3.37-4.87 

+ uncorrelated systematics 

2.84-7.84 

2.54-3.68 

3.55-9.69 

3.14-4.63 

Total Na reduction from 1st row 

33-36% 

24-35% 

39-40% 

30-39% 
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FIG. 6: As in Fig. 5, but for polynomial and uncorrelated systematic errors doubled (left) or halved (right), while statistical 
errors and systematics related to oscillation, normalization and resolution uncertainties are assumed to be the same. 


Figure 6 reports results analogous to the second half of Fig. 5, but with both the correlated polynomial and 
uncorrelated systematic uncertainties doubled, at the level of 3% (left) or halved, at the level of 0.75% (right). The 
left panels show that, with shape systematics at the few percent level, the hierarchy sensitivity tend to saturate in 
time, and can be lower than ~ 3tT in the worst cases, even after 10 years of data taking. Conversely, the right panels 
show that, with subpercent shape systematics, the sensitivity remains safely above 3cr (after 10 years) in all cases. 

Table II is analogous to Table I, but refers to the case where correlated polynomial and uncorrelated systematic 
uncertainties are doubled, as in Fig. 6 (left). By construction, the first two numerical rows are identical to Table I, 
while the others show a more pronounced reduction of the PINGU sensitivity, up to a factor of two when all the errors 
are included. 

Finally, Table III refers to the case where correlated polynomial and uncorrelated systematic uncertainties are 
halved, as in Fig. 6 (right). In this case, such uncertainties do not play a relevant role, and the reduction of the 
sensitivity (up to about 20-30%) is mainly due to resolution systematics. 

In conclusion, the results of Figs. 5 and 6 and of Tables I-III suggest that spectral shape uncertainties require a 
careful investigation, since they may be able to lower the PINGU sensitivity from 20% to 50%, as compared with an 
analysis including only the most obvious systematics due to oscillation and normalization uncertainties. 


TABLE II: As in Table I, but with correlated polynomial and uncorrelated systematic uncertainties taken at the level of 3%. 


5-year sensitivity No- 10-year sensitivity No- 


Errors included in the fit 

True NH 

True IH 

True NH 

True IH 

Stat. -1- syst (osc.-l-norm.) 

4.23-12.3 

3.34-5.64 

5.82-16.1 

4.49-7.64 

-I- resolution (scale, width) 

3.31-9.76 

2.95-4.37 

4.54-12.9 

4.00-5.94 

-I- polynomial (linear) 

3.03-8.79 

2.77-3.94 

4.07-11.5 

3.68-5.16 

-I- polynomial (quadratic) 

2.77-7.48 

2.46-3.60 

3.58-9.72 

3.13-4.72 

-I- polynomial (cubic) 

2.70-7.34 

2.41- 3.46 

3.43-9.32 

3.02-4.40 

-I- polynomial (quartic) 

2.66-7.27 

2.38-3.40 

3.36-9.21 

2.98-4.29 

-I- uncorrelated systematics 

2.37-6.48 

2.12-3.12 

2.75-7.50 

2.47-3.72 

Total Na reduction from 1st row 

33-53% 

37-45% 

46-47% 

45-49% 
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TABLE III: As in Table I, but with correlated polynomial and uncorrelated systematic uncertainties taken at the level of 0.75%. 


5-year sensitivity No- 10-year sensitivity No 


Errors included in the fit 

True NH 

True IH 

True NH 

True IH 

Stat. -I- syst (osc.-l-norm.) 

4.23-12.3 

3.34-5.64 

5.82-16.1 

4.49-7.64 

-I- resolution (scale, width) 

3.31-9.76 

2.95-4.37 

4.54-12.9 

4.00-5.94 

+ polynomial (linear) 

3.24-9.51 

2.92-4.30 

4.39-12.4 

3.92-5.77 

+ polynomial (quadratic) 

3.19-9.12 

2.84-4.16 

4.25-11.7 

3.74-5.47 

+ polynomial (cubic) 

2.18-9.07 

2.84- 4.14 

4.22-11.6 

3.73-5.43 

+ polynomial (quartic) 

3.15-8.94 

2.81-4.09 

4.17-11.5 

3.69-5.34 

-I- uncorrelated systematics 

3.11-8.86 

2.78-4.06 

4.08-11.3 

3.61-5.26 

Total No reduction from 1st row 

26-28% 

17-28% 

29-30% 

20-31% 


V. INTERPLAY BETWEEN HIERARCHY AND sin^ 6*23 DETERMINATIONS 

In the previous Section, we have assumed a prior range sin^ 023|true G [0.4, 0.6], roughly corresponding to the current 
±2cr allowed region [5]. In view of future constraints on 623 coming from ongoing and future accelerator experiments, 
it is useful to consider also a possible reduction of this range in prospective PINGU analyses. 

In particular, let us consider Fig. 7, which is analogous to Fig. 5, but is obtained for sin^ 6 > 23 |true G [0.46, 0.54]. 
One can notice a significant narrowing of the bands, and an overall gain in the minimum sensitivity. However, 
the pattern of progressive reduction of due to the inclusion of various shape systematics is similar to the one in 
Fig. 5. Therefore, prior information on sin^ 623 is of crucial relevance in determining the absolute sensitivity to the 
hierarchy, although it does not affect the relative reduction effects of systematic shape uncertainties. In this context, 
it makes sense to study how well this mixing angle may be determined by PINGU itself. 
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FIG. 7: As in Fig. 5, but for sin^ 023|true G [0.46, 0.54]. 
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FIG. 8 : Fitted value sin^ Sfs (at 1 , 2 and Scr) versus the true value sin^ ^ 23 ^'*, for the four possible cases where the test hierarchy 
(i.e., the one assumed in the fit) is either the true or the wrong one: (a) NH = true, NH = test; (b) NH = true, IH = test; (c) 
IH = true, IH = test; (d) IH = true, NH = test. 


Figure 8 shows, in each panel, the fitted value sin^ flfg (at I, 2 and Scr) as a function of the true value sin^ ^* 2 ™® G 
[0.4, 0.6], for the four possible cases where the true and tested hierarchies coincide or not. The results are obtained in 
a representative scenario with 5 years of PINGU data, and with polynomial and uncorrelated shape errors at the 1.5% 
level. In order to understand qualitatively such results, we stress that most of the hierarchy information (via matter 
effects) and octant-asymmetric information is embedded in the /i O e flavor oscillation channel, whose amplitude 
grows with sin^ 023- This information is enhanced in normal hierarchy, where matter effects are stronger for neutrinos, 
characterized by a larger cross section than antineutrinos. 

In Fig. 8 , the panel (a) refers to the case with true normal hierarchy, assumed to be unambiguously determined 
by PINGU. In this case, by construction, the fitted value sin^ flfg coincides with of the true value sin^ the 1 , 2 
and 3(7 bands provide then the accuracy of the sin^ 6*23 measurement with PINGU data only. The panel (c) shows 
analogous case for inverted hierarchy, which clearly shows a worsening of the accuracy for sin^ 023 j ^ much more 

pronounced effect of the octant degeneracy. The panel (b) refers to the case where the true hierarchy is normal, but 
PINGU is assumed to mis-identify it as inverted. In this case, the fitted value sin^ 0^3 is systematically higher than 
the true one; the reason is that the fitted IH spectrum tries to reproduce the intrinsically larger effects present in the 
true NH one, by increasing sin^ 0 I 3 as much as as possible. For analogous reasons, the opposite situation occurs in the 
panel (d), where the true hierarchy is inverted, but PINGU is assumed to mis-identify it as normal: the fitted value 
sin^ 0^3 is then systematically lower than the true one. Therefore, until the hierarchy is unambiguously determined, 
the determination of sin^ 023 may be subject to strong biases in PINGU. 

In conclusion. Fig. 7 illustrates the importance of prior information on sin^ 023 in determining the PINGU sensitivity 
to the hierarchy, while Fig. 8 illustrates, vice versa, the importance of prior information on the true hierarchy in 
determining the PINGU sensitivity to sin^ 023-^^ 


In this context, the role of the unknown phase <5 is marginal: we have verified that the reconstructed values of <5 are never constrained 
above the Icr level, for any of the PINGU error configurations examined in our fits (not shown). 
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VI. CONCLUSIONS AND PERSPECTIVES 

In this work we have examined, in the context of the proposed large-volume detector PINGU, several issues arising 
in the calculation of energy-angle distributions of atmospheric muon and electron events, and in the associated 
error estimates. In particular, it has been shown that the imprint of the neutrino mass hierarchy (either normal or 
inverted) on these spectra is sensitive to spectral shape variations at the level of (few) percent. This level of accuracy, 
not usually needed in fits to the dominant oscillation parameters (|Am^|,^ 23 ), poses unprecedented challenges to 
atmospheric neutrino experiments probing subdominant effects with very high statistics. 

In principle, one should try to characterize, at the percent level, all the known independent sources of uncertainties 
coming from models of differential atmospheric fluxes and cross sections, as well as from detector models used for event 
reconstruction. Breaking down such uncertainties into many separate nuisance parameters is crucial for a reliable 
estimate of the spectral “flexibility,” which unavoidably affects the hierarchy sensitivity in prospective or real data 
fits. Furthermore, in the limit of very high statistics, one should also account (via educated guesses) for residual, 
poorly known correlated and uncorrelated systematic uncertainties, which may escape a well-defined parametrization. 
Although each of these error sources may contribute to a tiny reduction of the hierarchy sensitivity, their cumulative 
effect may become noticeable. 

In this work, we have first analyzed the PINGU sensitivity to the hierarchy in the presence of the most obvious 
systematic errors, due to oscillation parameter and normalization uncertainties. Then we have added, in sequence, 
plausible shape systematics related to the resolution functions in energy and angle, generic polynomial shape deviations 
at the (few) percent level, and possible uncorrelated systematic errors at a comparable level. We have shown that 
their cumulative effect can induce a non negligible reduction of the PINGU sensitivity to the hierarchy—a result which 
deserves further studies, in order to reach more rehned and realistic error estimates. Finally, we have also discussed 
the interplay between the PINGU sensitivities to the hierarchy and to 023- 

In the context of atmospheric neutrino physics, a research program aiming at a better evaluation of the shape 
uncertainties of energy-angle spectra would be beneficial not only for PINGU, but for any future high-statistics 
atmospheric experiment, where such spectra may either embed subleading oscillation signals or provide a background 
for other kinds of emerging signals. In particular, all the issues discussed herein are expected to become more severe 
at lower (sub-GeV) neutrino energy thresholds, such as those proposed to probe GP violation effects |3D] . 

Rehned characterizations and evaluations of (known and unknown, correlated and uncorrelated) spectral uncertain¬ 
ties have been already considered in other physics contexts, including hts to parton distribution functions [22112,1] and 
precision cosmology data [5^ . They are now starting to be required also in the analysis of short-baseline reac¬ 
tor neutrino data [H]. We understand that the PINGU collaboration has recently initiated dedicated investigations 
towards similar objectives, with encouraging results [SUES. We think that, in the era of high-statistics precision 
experiments, these efforts should be largely promoted in the neutrino physics community, as it was recognized long 
ago pS] . 
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